ON THE CONVERGENCE OF LOCAL EXPANSIONS OF 
LAYER POTENTIALS 
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Abstract. In this paper, we study the convergence behavior of local expansions of singular or 
weakly singular layer potentials defined on a smooth, closed boundary T, when evaluated on T itself. 
For suitably chosen centers, we derive estimates for the rate of convergence of the induced local 
expansions as the centers approach the boundary. These estimates provide the analytical foundation 
for a recently developed quadrature method called 'quadrature by expansion' (QBX). They may also 
be of mathematical interest, particularly for asymptotic analysis in potential theory. 
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1. Introduction. The paper [16] describes a new method for the evaluation of 
layer potentials referred to as 'quadrature by expansion.' This method, denoted by 
QBX, is a technique for evaluating layer potentials on surface, which has simplified 
the development of fast, high-order accurate solvers for boundary integral equations. 
It is particularly easy to implement because it does not require the evaluation of 
integrals with non-smooth integrands. Here, we present the analytic foundations for 
QBX, which consist mainly in the establishment of decay estimates for one-sided 
local expansions induced by the layer potentials themselves. While introduced in 
the present context for the purpose of numerical analysis, these estimates may be of 
interest in their own right in asymptotic analysis and potential theory. 

For simplicity we assume that T C R n is a smooth compact hypersurface sepa- 
rating R" into two components R™ \ T = D + U D_. Throughout the paper, we let 
£>_ denote the bounded component and D + the unbounded component. By a layer 
potential, we mean an integral of the form 

F±(x) = ( k(x,y)ip(y)dS(y), for x e D±, (1.1) 



where ip(y) is a smooth density on T and k(x,y) is the Schwartz kernel, defined 
on R" x R™, for a pseudodifferential operator, IC, which satisfies the transmission 
condition, see Section 18.2 in [14]. We are primarily interested in the case where 
k(x,y) is the Green's function for the Laplace or Helmholtz equation (or the Cauchy 
kernel), but more general cases also arise in a number of applications. 

The kernels k{x,y) in these calculations have singularities on the diagonal (when 
x = y) which have hampered the development of stable, efficient, and high order 
accurate methods for their evaluation, particularly on surfaces in R 3 . The most fre- 
quently used high-order methods to date that are suitable for use with patch-based 
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surface discretizations have largely been based on product integration. In this ap- 
proach, the integral of the kernel multiplied by a piecewise polynomial approximation 
of the density <p on a piecewise smooth approximation of the boundary T is computed 
to high order using a mixture of analysis, linear algebraic techniques, and optimiza- 
tion. These have been developed in both two and three dimensions for a variety of 
kernels [6] [13j [25j [28]. Other methods, more closely related to QBX, are based on 
regularization of the kernel, combined with smooth quadrature rules and asymptotic 
correction [U [TU] HH HH1 [21] ■ Remarkably little use has been made, however, of the 
fact that the functions F± (x) have smooth extensions to D + and D_ . The limits from 
the two sides, of course, often differ. 

Remark 1. An exception is JP[ \15[ \20\j , in which the authors do make use of a 
local expansion about an off-surface point. In those papers, however, the viewpoint is 
global. In essence, a single expansion center is introduced (for Cauchy integrals with 
analytic data), with a radius of convergence determined by the location of the nearest 
singularity of the analytic function itself. 

In our case, we consider boundary data of finite differentiability, which is not the 
restriction of a real- or complex- analytic function and hence make no assumptions 
about the location of the nearest underlying singularity. We are interested in the 
behavior of local expansions in balls centered at points close to the boundary T, with 
radius equal to the distance from the center to the boundary. Due to the nature of 
the kernel functions, k(x,y), we can establish error estimates, valid in the closed ball, 
that only depend on the smoothness of the data and boundary near to the intersection 
of r with the boundary of the ball. 

A more detailed review of existing approaches to quadrature is given in [16j and 
the textbooks [3 03 [IS]. 




Fig. 1.1: A QBX expansion near a curve T, with a source point y, a target point Xq, 
and an expansion center x c . The angles <f> an d 9 are mainly used in Section [2731 

The QBX approach is based on the fact that, because of its one-sided smoothness, 
F±(x) supports a Taylor series expansion about x. More precisely, suppose that 
.To G r, that x c G £)_ and that there is a ball of radius r about x c , B r (x c ) C £)_, 
such that xq G bB r (x c ), where bB r (x c ) denotes the boundary of the ball. In other 
words, xo is a point of tangency of the ball and the surface T. (The analysis for 
the +- case is essentially identical.) The fact that F- is smooth in B r (x c ) and the 
remainder theorem for Taylor series shows that 

<M k (F_)r k+1 . (1.2) 
Here a denotes a standard multi-index in n dimensions, and Mk(F^) is a constant 



F_ (x - V -2 j^ Zc-ZO 

{a:\a\<k} 
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depending on the C fe+1 -norm of _F_ in B r (x c ). Because K, satisfies the transmission 
condition, this quantity is bounded by the Holder C fc+2+m ' Q -norm of if, where the m 
is the order of /C, and a is any positive number. The two key observations under- 
lying QBX are that the error in the expansion (|1.2[) can be controlled and that the 
derivatives {d"F_(x c )} can be computed accurately using only standard quadrature 
techniques for smooth functions. 

Briefly stated, QBX consists of a three-step procedure: For each boundary point 
xq, we 

1. find a nearby off-surface point x c , 

2. construct a local expansion of the form (11.2[) about the center x c of the held 
induced by the layer potential, and 

3. evaluate a partial sum of the local expansion from the previous step at xq- 
For kernels that satisfy a classical PDE, more efficient series expansion are preferable, 
based on standard separation of variables. For example, if K is a single or double 
layer potential for the Hclmholtz equation in three dimensions, for which the kernel is 
k(x, y) = e %k \ x ~y\ /{Att\x — y\), then it is convenient to use a (doubly-indexed) spherical 
harmonic expansion instead of a (triply-indexed) Taylor series. 

We restrict our attention in this paper to kernels k(x,y) which are Green's func- 
tions of some constant-coefficient elliptic partial differential operator, which are well- 
known to have their singular support along the diagonal. The analysis below, there- 
fore, applies to rcctihable hypersurfaces T, provided that T is smooth (C°°) in a 
neighborhood of the point, x, where we are trying to estimate the value of F±(x). 
To simplify the discussion below we assume that T is globally infinitely differentiable. 
This is not necessary; from our analysis it is clear that the accuracy of this approach 
near a point Xq on the boundary depends only on the smoothness of the boundary 
and the data near to xq. 

The purpose of this paper is to provide rigorous error estimates for the QBX 
procedure outlined above. In Sections|2]and|31 we derive decay estimates for analogues 
of Mk(F—) in p.2p for a selection of commonly encountered kernels. In Section dj we 
couple this analysis with quadrature error estimates for the numerical computation 
of {d£F-(x c )}. The corresponding analysis for the case of a general kernel is omitted. 
The steps involved are essentially the same as for the cases studied here. 

2. The two-dimensional case. We first consider two dimensional examples 
arising in potential and scattering theory, namely the Cauchy kernel, the Green's 
function for the Laplace equation, and the Green's function for the Helmholtz equa- 
tion. 

2.1. The Cauchy Kernel. Let D_ be a bounded domain in C = K 2 with a 
smooth boundary T — bD^, and ip a function on T. The Cauchy formula defines an 
analytic function / in C \ T : 



We denote the restrictions /_ = / [_d_ an d f+ = f \d + and recall the well-known 
jump formula: 

[/+ - /-] \r= f. (2.2) 
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The crucial fact in our analysis is that the functions f± extend to the closures of 
their domains of definition to be essentially as smooth as if. 

Remark 2. If smoothness is measured in terms of L 2 -Sobolev norms, then f± 
have a half- derivative more than, ip, whereas, if we use standard Holder norms, then 
f± have the same regularity, up to the boundary as the boundary data. Of course 
f± are analytic functions in the interior, so here we are primarily speaking of the 
boundary regularity. 

For purposes of definiteness, we suppose that z c E £>_, that the disk B r (z c ) C D-, 
and that zq E bB r (z c ) n T, where zq = z c + re l6 ° . In this disk: 

f-W = Y,—A V^iz-ZcY, (2.3) 
3=0 J ' 

where f^'(z) — d^f(z). In general this power series converges in B r (z c ), but in no 
larger disk. If tp has more than a half an L 2 -derivative, or is Holder continuous, then 
the series representation for /_ converges uniformly on B r (z c ). 

In the QBX approach, we would like to approximate the integral /_ from (|2.1|) . 
at the boundary point z$, by using a finite partial sum of the Taylor series (|2.3p . 
rather than a quadrature rule that is designed to handle the singular kernel head-on. 
Because of the special structure of the Cauchy kernel, two different approaches are 
available to estimate the error: 

e N (z) = /_(*) - —T 1 ^ ~ *c) j , (2.4) 
3=0 ] ' 

for z — z c + re l6 " . The first method does not generalize, but gives a simple explanation 
as to why this approach works, whereas the second gives a clear path to an error 
estimate in the general case. 

Let w(t) be an arc-length parametrization for T. The Cauchy integral is then 
given by 

1 ( , 5) 

Zm J w(t) — z 


To simplify notation, let us assume (without loss of generality) that z c = 0. For 
z E B r (Q), by considering the series for (1 — z/wft))" 1 , we obtain 

3 „,,// 



As 



then we can integrate by parts to obtain 

L 

/ x Z \ - 

e N 



W ' {t) --^0' (2.7) 



2m <rzLiJ V w W/ 3 -wit) 



j=N+l 
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In order to integrate by parts again, we use the fact that, for j > 1, 



-1 



-d. 



1 



[ w (t)Y (j-i)W(t) 'VN*)]'"- 1 

If we define the first order differential operator D t by 

g(t) 



D t g = d t 



w'{t) 



then we can integrate by parts N more times to obtain that 

L 



e N (z) 



2ni 



53 [ D?dt[<p(w(t))] 

— AT i n 



j=N+l 



z y- N -\j-N-iy. 



wty 

M<i 



w(t)j\ 



dt. 



(2.9) 



(2.10) 



(2.11) 



For N > 1, it is then straightforward to see that there is a constant, Mp ^, 
depending only on T and N such that 



(2.12) 



for any z with |z| < r. This proves 

Theorem 2.1. Let T be a smooth, bounded curve in C such that B r (0) C T c , but 
re lS ° E Tn bB r (0), where T c denotes the complement of T C R 2 . For eac/i positive 
integer N there is a constant Afr.w suc/i £/iai, /or 6 C w+1 (r), £/ie error in the 
truncated Taylor series approximation is given by 



1 f ip(w)dw 
2ni I w — re l9 ° 



2^— f -7, e ' 



< Mr : Ar||(p|| C Jv+i (r )r 



JV+1 



(2.13) 



Note that the coefficients of the expansion can be obtained by evaluating the non- 
singular integrals: 



/«(0) 



1 



ip(w{t))w' (t)dt 



j\ 2m J [w{t)Y +1 



(2.14) 



Error estimates for the numerical computation of these integrals are discussed in 
section H) 

We turn now to a second approach for estimating the error in truncating the 
Taylor series, which uses the Fourier expansion of /_ on bB r (Q) directly. As /_ is 
holomorphic in J5 r (0) its Fourier series on the boundary only has terms with non- 
negative exponents: 



3=0 



where 



2tt 



(2.15) 



(2.16) 
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We again want to show that 



/_(re")-5>. 



•jo 



3=0 



0{ 



(2.17) 



which is evidently a matter of estimating the Fourier coefficients themselves. 
Integrating by parts, and assuming j > 0, we have 



27T 



(2.18) 



Note that e has been replaced by e *w which is a consequence of the fact 
that 



dof-(re t8 ) = ire w f W (re w ). 



(2.19) 



This shift in degree reappears when we consider spherical harmonic expansions below. 

If /_ has N + 1 derivatives on bB r (0), then we can repeat this (N + l)-times to 
obtain that, for j > N + 1, 



1 

2ni 



JV+l 



(j-N- 1) 



2tt 



f {N+1 \re w )e-^- N ^ e de, 



(2.20) 



from which the estimate in (|2.17l) is immediate. The implied constant depends on 
||/_||cJv+i(s r (o))) which we know from the mapping properties of the Cauchy integral 
operator is of the order 0(||y > ||c JV + 1 ^(r)) f° r an Y P > 0, see Theorem 1.4 on page 147 
of [27] ■ We therefore get an estimate slightly weaker than that obtained in (|2.13l) . 

2.2. Harmonic Layer Potentials. We now apply similar ideas to the evalua- 
tion of single and double layer potentials on the boundaries of planar regions. The 
fundamental solution for the Laplace operator is given by 



G (z,w) = ^-log|z 

Z7T 



— Relog^-to). 

Z7T 



(2.21) 



While the complex log is multi-valued, its real part and its complex derivative are 
globally defined as single- valued analytic functions. As before we let D_ be a bounded 
region with boundary T. Assuming that ip is real valued, the single layer potential 
with density supported on T is defined as 



1 

2^ 



Re 



<p(w{t))log(w(t) - z)dt 



(2.22) 



where z = + and w(t) is again an arc-length parametrization of T, viewed as 
lying in the complex plane. 

We assume (without loss of generality) that the origin E D_ , which we take to 
be the center of our power series expansion. If B r (0) C D_, then for z £ B r (0) we 
can rewrite this integral as 



Z7T 



<p(w(t)) log 



1 - 



w(t) 



dt 



A , 



(2.23) 
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where 



A 



2n 



¥>(«>(*)) log \w(t)\dt. 



(2.24) 



Let us denote by u\ the first term in (|2.23p . As \z/w(t)\ < 1, the power series 
expansion for the principal branch of log(l — £) about C = 0, shows that u\ has the 
following power series expansion: 



i 



1 ( z 



;(t) 



(2.25) 



To approximate it(£i,£2) for z € -Br(O) we use Ao plus a finite partial sum of this 
series, giving an error term: 



ejv(£i)6i) = u(Ci,6) 



I ^ JV 



«;(<) 



1 f z 



j \w{t) 



dt 



^ ]T Rc / r -(/ri/i)- 
Using the identity in equation (|2.9|) we can integrate by parts A-times to obtain: 

L 



(2.26) 



z N+1 (j-N)\ 
ejvlti,t2) = — > ^ -Re 



.7 



w(t) 



j-N-1 



j=N+l 

where D t is again the differential operator defined in (|2 . 10|) . As r —> 0, 

dt 



dt 



oc — log r. 



(2.27) 
(2.28) 



The preceding estimates yield the following theorem. 

Theorem 2.2. Suppose that T is a smooth, bounded curve embedded in R 2 , such 
that B r (0) C r c . but zq — re 10 " G bB r (Q) n T. For N a positive integer, there is a 
constant Mjv,r so that if ip G C N (T), then 



Go(zo,y)f(y)ds(y) - 



Ao-^ReU^))tj(^ydt 



<M NX M CN{r) r N+1 \og-. (2.29) 



The constant Aq is defined in (|2.24l) . 

The double layer potential is defined as 



"(£1,62) = / d nw log\w - z\ip(w(t))dt, 



(2.30) 
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where n is the outward normal along 6_D_. As log(io — z) is analytic, the Cauchy- 
Riemann equations imply that 

dn m log \w — z\ = Re[dn w \og(w — z)] = — lm[d Tw \og(w — z)], (2.31) 

where r is the unit tangent vector along TD_. Thus, we may write 

d nw log\w-z\ = -Im ■ (2.32) 

w{t) — z 

It follows that the estimate proved in the previous section for the Cauchy transform 
can be applied to analyze the error when replacing the double layer potential, , £°) 
with the finite partial sum: 

1 N 

n j=o 

where (r, 9 ) are the polar coordinates of (£i,£;>) with respect to the origin (which is 
the center of the local expansion) . 

Corollary 2.3. Let T be a smooth, bounded curve in C such that B r (0) C T c , 
with zq = £i+*£2 = re 10 " £ TP\bB r (0). For each positive integer N there is a constant 
Mjv,r such that if ip G C N+1 (T), then 

HZi,$)-v {N )(r,e )\ < M Nir r N+1 \\<p\\ c » +Hr) , (2.34) 

where v is defined in (|2.30|) and vr m is defined in ()2.33j) . 

2.3. Layer Potentials for the Helmholtz Equation. Our last example in 
two dimensions concerns layer potentials that satisfy the Helmholtz equation: 

(A + fc 2 )it = 0. (2.35) 

The fundamental solution, for k ^ 0, is given by the zeroth order Hankel function of 
the first kind: 

G k {x,y)= % -Hi 1 \k\x-y\). (2.36) 

It is well-known that the kernel Gk(x, y) defines a classical pseudodifferential operator 
of order —2, see [26]. For our purposes, the most important fact about this function 
is the Graf addition theorem, see equation 9.1.79 in pQ, which gives the series repre- 
sentation: 

00 

H^(k\x-y\)= J\i\(k\x\)H^{k\y\)e^ B -*\ (2.37) 
i— — 00 

valid so long as 

\x\ < \y\ where x = \x\e t6 and y = |y|e 10 . (2.38) 

Here Ji and denote the usual Bessel and Hankel functions of the first kind, 

respectively. 



re ie °y ip(w(t))w' (t)dt 



,{t) 



w(t) 



(2.33) 
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We assume (without loss of generality) that the origin E D_ , which we take to 
be the center of our series expansion. With these assumptions about D_ and B r (0), 
we define the single layer potential: 



u ( x ) = J <p(y)Gk(x,y)ds(y) 
r 

= £ J W (k\x\)e iW / H^(k\y\)e- a ^ V {y)d S (y) 



(2.39) 



a i J \i\( k \ 



x e 



r 

no 



l=-c 



The second equation follows from (|2.37[) The series in the last line is valid for x £ 
B r (0), with the coefficients {a;} defined by the integrals: 



ai = I Hff(k\y\)e- U *»<p(y)ds(y) 



l<l>y 



(2.40) 



We suppose that tp e C Ar " 3 (r), for some (3 > 0, so that u± belongs to C N+1 ^{D±). 
If we assume that xo = re l6 " £ bB r (0) fl T, then, in the QBX approach, we would like 
to use the partial sum of the series: 



N 



^ aiJ\i\(kr)e 



iW 



(2.41) 



l=-N 



as an approximation for u-(xq). Because the Fourier series representation is unique, 
it follows that 



2tt 



otiJ\i\(kr) — — I u- (r cos 9, r sin 9)e 



-no 



(2.42) 



We now seek, as above, to estimate the error 



TV 



(-(xo)- ]T «,J|,|(Ar)e <Wo 



l=-N 



(2.43) 



We make use of the integration by parts formula: 

2tt 

— / U-(r cos 9, r sin 9)e~ 



-H8 



2tt ■ I 



where 



[dgU-] (r cos 9, r sin 9)e 



-i(/+i)e 



[d z uJ\(rcos9,rsm9)e- l{l - 1)e ^ d9, (2.44) 



d z = ^(d x + idy) and d g = -(d x - id y ). 



(2.45) 
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If we integrate by parts, in this manner, N + 1 times, then the order / term, where 
./V < |Z|, produces an integrand of the form: 



N+l 



1 J2 ci,N+i,je^ l+2j - N -^ 6 [did^ +1 - j u^}(rcos9,rsiae). (2.46) 



3=0 

Each coefficient ci t N+x,j is a sum of 



N + l 
j 



terms of the form: 



1(1 + ei) • • • (I + e N ) 
with each satisfying the upper bound: 

R,iV+l,ji < 

I'l' 

Thus, for \l\ > N, summing over j gives the estimate 



where Cj G {-N, 1 - N, . . . , N - 1, N}, (2.47) 



\l\-N-l)\ 



otiJu\(kr)\ < 



(2T-) iV + 1 (|Z| — TV — 1)! • |h 



Hence, for N > 1, there is a constant Mn so that 

N 



L-{re 100 )- <xiJ\i\(kr)e iWo 



l=-N 



<M N r N+1 \\u-\\ cN+1 



e N + 1 (B r {o)y 



(2.48) 



(2.49) 



(2.50) 



The operator G k maps C N ^(T) to C N+ ^(DJ) for any /3 > 0, see [TJ. Hence the 
preceding estimates yield the following result. 

Theorem 2.4. Suppose that T is a smooth, bounded curve embedded in R 2 , such 
that B r (0) C T c , with re 100 G bB r (0) n T. For k G [0, oo), N a positive integer, and 
(3 > 0, there is a constant M' N ^ so that if ip G C JV, ' 3 (r), i/ien 



/ G k (re ie °,yMy)ds(y)- ]T a,J|,|(fr 

" / jvr 



r)e 



l=-N 



< 



M' Nt pr"+ l \\<p\\ cN , Hr) . (2.51) 



i7ere i/ie coefficients {a;} are given by (|2.40[) . 
An estimate for the double layer 



v(x) = J d ny G k (x,y)p(y)ds(y) 
r 



(2.52) 



is obtained just as easily. In £)_, the double layer has an expansion of the form 

oo 

«_(*)= E a^|z|(^M)e^, (2.53) 



Z = — OO 
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valid for x £ £> r (0), where now the coefficients {a{\ are defined by the integrals: 

m = I ' d ny [H^{k\y\)e' a ^My)ds(y). (2.54) 



Arguing exactly as before we can show that, for N > 2, we have the estimate: 

N-l 

IS ^- 71/f „N\ 

~M<J\l\\Kyj;\)ti 

1=1- N 



v(x Q ) - ^ anJ\i\(k\x\)e' 



< M N r N \\v-\\ CN(s - r 



(o))- 



(2.55) 



Because the double layer defines an operator of order —1 we have 

Theorem 2.5. Suppose that T is a smooth, bounded curve embedded in R 2 , such 
that B r (0) C T c , but re m € bB r (0) n T. For k £ [0, oo), N a positive integer, and 
{3 > 0, there is a constant Mpf a so that if ip G C N ' l3 (T), then 



r N-l 

d ny G k (re ie °,y)p{y)ds(y)- £ a, J\i\(kr)e iW ° 

J 1 1 AT 



1=1- N 



< 



M'^McN.ep)- (2.56) 



Here the coefficients {cti} are given by (|2.54l) . 



3. Three-dimensional Layer Potentials. In this section we consider the eval- 
uation of layer potentials arising from the Helmholtz equation in three dimensions. 
The "outgoing" fundamental solution is given by 



Gk(x,y) = 



e ik\x-y\ 

An\x — y\ 



(3-1) 



for Irak > 0. The relevant formulae in the harmonic (fc = 0) case are somewhat 
different, but follow the same lines as those for k ^ 0. We leave that analysis to the 
interested reader. 

As in the two-dimensional case, this Green's function also has a classical expansion 
in terms of spherical harmonics [21 : 



G k (x,y)=ikJ2j l (k\x\)h l (k\y\) ]T Y™ (6\ , <j)i)Yf m (# 2 , <ht) 



1=0 



m= — l 



assuming |x| < \y\, with 



where 



x = \x\u(9i,<j>i) and y = \y\us (9 2,^2) 



(sin (j> cos 9, sin (j> sin 9, cos < 



(3.2) 



(3.3) 



(3.4) 



Here, ji and hi denote the spherical Bessel and Hankel functions of the first kind, 
respectively [1 . The normalized spherical harmonics YJ m are defined for |m| < I by 



17 



4tt (Z + |mQ! ' v ' 



(3.5) 
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where P™ is the associated Legendre function of degree I and order m. 

We let D_ denote a bounded region in M 3 with smooth boundary T. We assume, 
as in the two-dimensional case, that the origin G D_ and that B r (Q) C with 
x e 6B r (0)nr. If ip is an integrable function defined on T, then a single layer potential 
with this density is given by 



f(x)= / G k (x,y)<p(y)dS(y) 



(3.6) 



1=0 m=-l 

The series representation is valid in B r (0), with 



ai m =ikj h l (k\y\)ip(y)Y l - m (9 yi ^ y )dS(y). (3.7) 
r 

Once again, we would like to approximate f(xo) by the partial sum of this expansion: 

N I 

/(so) « / (JV) (^o) = X)j'i(*kol) E « im rr(^o>o). (3.8) 



;=o m=-l 



To estimate the error \f(xo) — f( N ' ) (xo)\, we use an argument very much like that used 
to obtain (|2.50[) and (|2.5ip . Before we turn to that task, however, let us compute the 
series expansion of the double layer 



<x) = J d ny G k (x,y)ip(y)dS(y). (3.9) 
r 

v(x) also has an expansion like that in the second line of (|3.6[) . 

oo / 

v{x) = X)ji(*N) Yl aimYr(Ou<f>i), for \x\ < r, (3.10) 

1=0 m=-l 

where the coefficients are now defined by the integrals: 

a lm =ik [ d ny [fi^yDY-™^,^)] v(y)dS(y). (3.11) 



The value of v(xo) can again be approximated by taking a partial sum of the series 
in (|3T0)l . 

The error estimates for both of these series approximations come from determining 
the r dependence of the coefficients in these expansion. We carry out the analysis for 
the single layer and leave the details for the double layer to the interested reader. 

The necessary estimates follow from integration by parts using the fact that: 



L±1T = V(iTm)(I±m + l)yi ra±1 , (3.12) 
where the vector fields are defined in (8, 0)-coordinates by: 

L± = e ±w [<9 ± i cot 4>d ] . (3.13) 
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An elementary calculation shows that, as operators on L 2 (S 2 ), the adjoints of L± are 
— L±. In integrating by parts, we need to apply these vector fields to g(ruj(0, </>)), 
where g is a function smooth in the closure of the sphere of radius r. An elementary 
calculation gives: 

L ± g{rio{9 ,(!))) = rcos^[ff x ± i 5y ](rw(M)) - re ±4e sin0 5z (™(#, 0)). (3.14) 
It is useful to observe that: 

cos0 = \ l^-Yi and e ±lf> sin</> = J^-Y^ 1 , (3.15) 

V o \ o 



and therefore 



L ± g(roj(9,<f>)) = rJ — [Y°(g x ± ig y ) - v^lf jj(ru(9,^)). (3.16) 



The Clebsch-Gordan relations [23] (or [22], 134.3.201 134.2.11 134.2.30 , show that 



there are collections of coefficients {ct . n ; c 7m . n } so that for j e {—1, 0, 1} we have 



the identities 



l-i i+i 

n^- m =Ew^-i+ E c k« y w- ( 3 - 17 ) 

n=l-i n=-(X+0 

An elementary computation shows that 

ii?'ir>)i 2 < 47r (i 3 +|j|) i y '" m Hi 2 - (3-i8) 

Integrating this inequality over the sphere, and using the orthogonality relations sat- 
isfied by the functions {Y™}, we easily deduce that for j E { — 1, 0, 1}, 

l-i l+i 3 

E IW| 2 + E l4m;J 2 < 47r(1 + | . |r (3-19) 

n=l-l n=-(2+l) V lJU 

As before, we observe that the representation f(roj(9, <f>)) of the single layer poten- 
tial in (|3.6p is unique, and therefore the coefficients can be computed by integration 
on a sphere of radius r about the origin: 

27T 7T 

ai m (r)=ji(kr)a lm =[ J f {ru{6 , ^Yf™ (9, 0) sin <t>d<j>de. (3.20) 
o o 

For simplicity assume that m > 0, and use (|3.12|) to see that 



L+Y- m = v / {l + m)(l-m + l)Y l 1 - m . (3.21) 
Integrating by parts in the integral defining S; m (r) we see that 

r r Y n (9 6) 

a lm (r) = - L + [f(ruj(9, <p))}—==J====J==dS{9, 4>) (3.22) 
J J + m){l -m + 1) 
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Using the relations in (|3 . 1 6[) and (|3.17[) . we see that the integrand becomes 




in {( f x + if y )Y? - V^fzY^ 1 - 
+ m)(l - m + 1) 



47T I 

3(l + m)(l-m+l) { ifx + lfv) 



E C 0lm-,rXi n -l + E C tlm;rXl+l 



n=X—l 



n=-(l+0 



l-l 



l+l 



E c llm;nXl^l + E C lZm;n^2 + l 
n=-(l+Q 



n=l-2 



(3.23) 



The formula in (|3.23|l has several notable features: 

1. The entire quantity is multiplied by r. 

2. The terms in these sums are of exactly the same types as those that appear 
in the original definition, with the integrals computing spherical harmonic 
expansion coefficients of derivatives of / along |x| = r. 

3. Starting with a spherical harmonic of degree I we obtain terms of degree I — 1 
and l + l, showing that this procedure can be repeated I times. 

4. The coefficients {cf lm . n } can be bounded by the identity in (|3.19p . 

Using these observations, we carry out integration by parts N + 1 times to obtain an 
estimate for the remainder: 



eW(*o) = /(*o)-/W(z )= E 



i 

E 

N+l m=-\ 



5b„ coir (*o, 6>). 



(3.24) 



Integrating the right hand side of (|3.23p over the unit sphere leads to the com- 
putation of certain spherical harmonic coefficients of (f x + if y ), and v2/ z restricted 
to the sphere of radius r. In order to estimate the coefficients {az m (r)}, we need to 
introduce some notation. We let 



Di=d x + id y , D 2 = d x - id v , D 3 = V2d z 



(3.25) 



For a 3-multi-index 7 = (71,72,73), we let {cEiZ, (r)} denote the spherical harmonic 
coefficients of D 1 f f^m), where D 1 is the constant coefficient differential operator 
DfDfDf, so that 



42 00 = 



[Z) 7 /](rw(6», </>))F, m (6», (ft) smcj)d(j)dd. 



(3.26) 



To obtain estimates for {(Xi m (r)} that are proportional to r N+1 , we need to inte- 
grate by parts in (|3.20j) N + l times. Each integration by parts increases the number 
of terms by a factor of 4. Whether one chooses L + or L_ at each step is not too 
important, but for the fact that for any pair (l,m) one can always choose the sign so 
that coefficient on the right hand side of (13.121) does not vanish. 

Our formula expresses 



a lm (r)Yr(0,cb) 



(3.27) 
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as a sum of (21 + 1) x 4 N+1 sums, each consisting of terms of the form 

r N+1 -a^ vi (r)C 1 ---C N+1 Y l u ^(9,4). 



47r\ 2 



Each matrix C q is of the form 



(3.28) 



S~i jqlqUq ',V q 

°? - J • 



(3.29) 



For each q, the pair u q , v q are the matrix indices in C, 



j q €{-1,0,1}, 
I - (N + 1) <l q < I + (N + 1), and therefore 
u q E {-l q , ...,lg} and v q E {-(l q ±l),...,(l q ± 1)}. 



(3.30) 
(3.31) 
(3.32) 



as 



We think of afcl (r) as a 1 x [2 (I i ± 1) + 1] vector (indexed by vi), and Y t 
a (2/ + 1) X 1 vector of functions (indexed by upr+i)- Moreover 7 is a 3-multi-index 
satisfying \ j\ = N + 1. 

Finally, the denominator is of the form 



d q = V ( l Q T u q)( l q ± Uq + 1), 



(3.33) 



with + or — chosen so that d q never vanishes. With this understood, d q is easily seen 
to satisfy the estimate: 



'2l q < d q . 



(3.34) 



8) 



Using p.!9p and (|3.34[) . we can bound the Frobenius norm of each matrix C, 
in (|3.28j) (with u q and v q as the matrix indices) by y / 3/47r. Hence the absolute value 
of the quantity in (|3.28j) is bounded by 



JV+1||C;(7) 



where we use the fact that pointwise 

l 



E i y r(M)i s 



-/ 



21 + 1 

4tt : 



21 + 1 
Air 



(3.35) 



(3.36) 



It therefore follows that for some pairs {(jj,Pj) ■ j = 1, (2Z + 1) •4 JV+1 } where 
= N + 1 and pj E { — (N + 1), . . . , (N + 1)}, chosen as described above, we have 



E *i™(r)Y? 



m=—l 



. (2/+l).4" + 1 



4?r 



(3.37) 



Finally, collecting terms we see that there is a constant so that: 

i 

i 



\eW( X0 )\<K N r^ E E 

76J"3,w+i (=0 



E ia£ } (0l a 



.m=—l 



(l + N + 1)^ 



(3.38) 
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Here J3 jv+i is the collection of 3-multi-indices with |-y| = TV + 1. In order for these 
infinite sums to be finite, we need to assume that, for each 7 G J3,n+i the sum 
satisfies: 



E ; l^(r)| 2 =o(^), (3.39) 



for some 6 > 0. 

The estimate in (|3.39p holds if, for each 7 g J3,n+i, the function D 1 f \bB r (o) 
belongs to the L 2 -Sobolev space W 3+5 ' 2 (bB r (0)). For a single layer, this would follow 
from assuming that ip G W 3+N+S ' 2 (T). This represents a small loss of regularity over 
the estimates we obtained in two dimensions. This loss would appear to be a result of 
the added complexity of multiplying spherical harmonics rather than exponentials. A 
more careful estimate of the terms in p. 281) might provide a somewhat better result. 

Theorem 3.1. Let k be a complex number with Imfc > and let T denote a 
smooth hyper surface in R 3 , such that B r (0) C T c with xq = ru>(0o, (f>o) € T n bB r (Q). 
For each positive integer N and S > 0, there is a constant M^,s such that 



r 



N I 

G k {x ,y)^{y)dS{y)~Y.^ kr ) E «^ m (^o,0o) 

1=0 m=-l 



< 



M N r N+1 \\cp\\ W 3 + N +l ,2 ir) , (3.40) 
so long as ip £ W 3+N+5 ' 2 (T). The coefficients {ai m } are given by (|3.7[) . 

4. Quadrature error. The prior estimates establish bounds on the tails of the 
expansions used by QBX, and hence on the error incurred in their truncation. They do 
not take into account that the coefficients in these expansions must still be computed 
numerically. In this section, we derive an error estimate for this quadrature problem. 
This estimate then helps to establish an upper bound on the "mesh spacing" h needed, 
leading to straightforward control of the overall error in QBX, which is bounded above 
by the sum of these truncation and quadrature errors. 

Concentrating for the moment on the single layer for the Laplace equation in two 
dimensions, let us first rewrite (|2.26[) as 

N 

e N (ti,&) =u(^,^)-A -Re{ YV-^ } =Re{ > ] V (4.1) 



j=l I I j=N+l 



where 



L 

A " = hJ <P(w(t))logHt)\dt, 


1 tvW)) 



2ir J j[w(t)Y 



i) 



and T = w([0,L)). Recall, also, that the center z c here is assumed to be the origin, 
for simplicity of notation. In practice, what we wish to estimate is actually 



ew(£i,60 



«(&,6)-4?-Re 



N 



(4.2) 
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where A® and af are the computed approximants of A and a,-. We write this error 



e at (£1,6) 



N 



- A - Re < 



A" 



- (A Q - J 4 )-Re<|^[ S Q -a J ]; 



(4.3) 



The hrst term is 0(r + ) as shown in Theorem 12.21 The second term depends 
strongly on the specific quadrature rule used. To facilitate straightforward generaliza- 
tion to complicated geometries in more than two dimensions, we choose a composite 
Gauss quadrature for this role. If the boundary T were divided into M equal subinter- 
vals of size h, and a gth order Gauss quadrature rule were used on each subinterval, 
then 



A? 



CJ\.T)( AV| M | cas . (4.4) 



To see this, we note that on a curve segment I\ of length h, the standard estimate 
for g-point Gauss-Legendre quadrature [8; yields 



y>(w(* ra ) ) 

li[^(tn)] 



< 



( g !) 4 

j(2g + l)(2«)!3 



(4.5) 



where D 9 denotes the q th derivative of the integrand with respect to the integration 
parameter along I\, and {t n ,w n } denote the Gauss nodes and weights scaled to IV 
From Stirling's approximation 



27m™+2 e -" < n ! < 2^Fn n+ V™, 



we have 



ip(w(tn)) 



< 



h 2 1 +1 

(2g)!4 2 « 



£> 2<? 



<?Mt)) 



'(*)]• 



Omitting a detailed computation of the derivative, we note that it satisfies a bound 
of the general form 



D 2 « 



<p(w(t)) 



[w(t)\ 



<c q (N,r)^\\<p\y q . 



The constant in this estimate depends on the order of the local expansion iV and on 
the smoothness of the curve T, through a polynomial in (w, u>W , . . . , w^ 2q%> ). Since the 
dependence of the constant on N, and the parametrization are rather complicated, 
but independent of ip, we leave it in the form above. 
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Adding up the errors from all subintervals yields (|4.4[) . The total error incurred 
by the QBX method can therefore be written in the form 

E = 0(r N+1 ) + C q (N,T) (AV^IIc*- 
Fixing the order of the expansion N and letting r — Ah yields an error of the form: 
E = 0(h N+1 ) + C q (N,T) flV||p|| cas . 

Implemented in this manner, the QBX method is asymptotic, converging like a method 
of order N+l, until the error is dominated by the second term. To obtain a convergent 
scheme one can, for example, set r = y/h so that 

E = 0{h^) + C q (N,T)(j£\ Hie*. 

Similar estimates can be obtained for Helmholtz layer potentials and for three-dimensiona 
cases. 

In practice, when dealing with complicated boundaries, adaptive discretization 
is generally required, and the assumption that the boundary is divided into equal 
subintervals must be relaxed. This does not change the error analysis in a substantial 
way, and we refer the reader to |16) for details and examples. 

5. Conclusions. It is well-known that local expansions of layer potentials are 
analytic away from the boundary. There are surprisingly few results in the literature, 
however, concerning the behavior of such expansions when evaluated at the limit of 
their radii of convergence - namely at the nearest boundary point itself. The estimates 
in Theorems 12.112.21 \2A\ I2.5[ and 13.11 show that, despite the absence of a geometric 
decay parameter, high order accuracy can be obtained in computing layer potentials 
in a manner that is controlled by the smoothness of the boundary data. The analysis 
presented here serves as a complement to |16j , providing a convergence theory for the 
quadrature method QBX. 

We note that sharper results can be obtained [3j , if some assumptions are made 
about analyticity of the data and the location of the nearest singularity. 

Finally, estimates of the type obtained here may be of interest in areas of mathe- 
matics that involve smooth continuation, microfocal analysis, and potential theory. 

Acknowledgment. The authors would like to thank Alex Barnett, Zydrunas 
Gimbutas and Michael O'Neil for many useful discussions. 

REFERENCES 

[1] M. Abramowitz and I. Stegun (editors), "Handbook of Mathematical Functions", National 

Bureau of Standards, 1964. 
[2] K. E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge 

University Press, Cambridge, UK, 1997. 
[3] A. Barnett, "Evaluation of layer potentials close to the boundary for Laplace and Helmholtz 

problems on analytic planar domains", submitted, 2011. 
[4] J. T. Beale AND M.-C. Lai, "A method for computing nearly singular integrals", SIAM Journal 

on Scientific Computing, 38, 1902-1925, 2001, |doi: 10.1137/ S0036142999362845| 
[5] C. A. Brebbia, J. C. F. Telles, and L. C. WrOBEL, Boundary Element Techniques, Springer, 

New York, 1984. 



ON THE CONVERGENCE OF LOCAL EXPANSIONS OF LAYER POTENTIALS 



19 



[6] J. Bremer, Z. Gimbutas, and V. Rokhlin, "A nonlinear optimization procedure for 
generalized Gaussian quadratures", SIAM J. Sci. Comput., 32, 1761-1788, 2010, 
[doi: 10.1137/080737046| 

[7] David Colton and Rainer. KR.ESS Integral equation methods in scattering theory, Krieger Pub- 
lishing Co., Malabar, Florida, reprint ed., 1992. 

[8] P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, San Diego, 
1984. 

[9] J. N. Lyness and L. M. Delves, "A numerical method for locating the zeros of an analytic 

function", Mathematics of Computation, 21, 543-560, 1967, |doi: 10.2307/200 4999 
[10] J. Goodman, T. Y. Hou, and J. Lowengrub, "High-order and efficient methods for the vortic- 

ity formulation of the Euler equations" , Communications on Pure and Applied Mathematics, 

43, 415-430, 1990, |doi: 10.1137/0914067| 
[11] L. Greengard, A. Klockner, and M. O'Neil, "Fast algorithms for the evaluation of layer 

potentials using 'Quadrature By Expansion' " Technical report, Courant Insitute, 2012 (in 

preparation). 

[12] D. J. Haroldsen and D. I. Meiron, "Numerical Calculation of Three-dimensional Interfacial 

Potential Flows using the Point Vortex Method" , Communications on Pure and Applied 

Mathematics, 43, 415-430, 1990. |doi: 10.1137/S10648275 96302060 
[13] J. Helsing and R. Ojala, "Corner singularities for elliptic problems: integral equations, graded 

meshes, quadrature, and compressed inverse preconditioning", J. Comput. Phys., 227, 

8820-8840, 2008, |doi: 10.1016/j.jcp.2008.06.022| 
[14] Lars Hormander, The Analysis of Linear Partial Differential Operators III, Springer- Verlag, 

Berlin, 1985. 

[15] N. I. Ioakimidis, K. E. Papadakis, and E. A. Perdios, "Numerical evaluation of analytic 
functions by Cauchy's theorem", BIT, 31, 276-285, 1991, |doi: 10.1007/BF01931287| 

[16] A. Kloeckner, A. Barnett, L. Greengard, and M. O'Neil, Quadrature by Expansion: A 
New Method for the Evaluation of Layer Potentials, arXiv 1207.4461 2012. 

[17] R. Kress, Numerical Analysis, GTM, Springer- Verlag, New York, 1998. 

[18] R. Kress, Linear Integral Equations, Applied Mathematical Sciences, vol. 82, Springer, 1999. 
[19] J. Lowengrub, M. Shelley, and B. Merriman, "High-order and efficient methods for the 

vorticity formulation of the Euler equations" , SIAM Journal on Scientific Computing, 14, 

1107-1142, 1993, |doi: 10.1137/0914067| 
[20] J. N. Lyness and L. M. Delves, "On numerical contour integration round a closed contour", 

Mathematics of Computation, 21, 561-577, 1967, |doi: 10.1090/S0025-5718-1967-0229388-0"1 
[21] MORSE and Feshbach, Methods of Theoretical Physics, McGraw-Hill, New York, 1953. 
[22] F. W. J. Olver (ed.), NIST Handbook of Mathematical Functions, Cambridge University Press, 

2010. 

[23] M. E. Rose, Elementary Theory of Angular Momentum, Dover, New York, 1995. 

[24] C. Schwab and W. L. Wendland, "On numerical cubatures of singular surface in- 
tegrals in boundary element methods", Numerische Mathematik, 62, 342-369, 1992, 
[doi: 10.1007/BF0139 6234 

[25] J. Strain, "Locally-corrected multidimensional quadrature rules for singular functions", SIAM 
Journal on Scientific Computing, 16, 992-1017, 1995, |doi: 10.1137/0916058| 

[26] Michael J. Taylor, Partial Differential Equations II, in series Applied Mathematical Sciences, 
vol. 116, Srpinger- Verlag, New York, 1996. 

[27] GuO-Chen Wen, Conformal Mapping and Boundary Value Problems in series Translations of 
Mathematical Monographs vol. 106, American Mathematical Society, Providence, 1992. 

[28] N. Yarvin and V. Rokhlin, "Generalized Gaussian Quadratures and Singular Value Decompo- 
sitions of Integral Operators", SIAM Journal on Scientific Computing, 20, 699-718, 1998, 
[doi: 10.1137/S1064827596310779| 



